#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <stdlib.h>

#include "raylib.h"
#include "pub_lib.h"
#include "FreeImage.h"
#include "image.h"
#include "dpi_algo.h"

/*
 * 图像反转
 * 图像数据来源于 ori, 修改后的数据写入 chg
 * */
void img_negative(FIBITMAP* ori, FIBITMAP* chg)
{
	unsigned int x, y;
	BYTE value;

	unsigned int w, h;
	w = FreeImage_GetWidth(ori);
	h = FreeImage_GetHeight(ori);

	for (y = 0; y < h; y++) {
		for (x = 0; x < w; x++) {
			FreeImage_GetPixelIndex(ori, x, y, &value);
			value = 255 - value;
			FreeImage_SetPixelIndex(chg, x, y, &value);
		}
	}
}

/*
 * 两级图像
 * 灰度级 < threshold, value = 0
 * 灰度级 >= threshold, value = 255
 * */
void image_binary(FIBITMAP* ori, FIBITMAP* chg, unsigned char threshold)
{
	unsigned int x, y;
	BYTE value;

	unsigned int w, h;
	w = FreeImage_GetWidth(ori);
	h = FreeImage_GetHeight(ori);

	for (y = 0; y < h; y++) {
		for (x = 0; x < w; x++) {
			FreeImage_GetPixelIndex(ori, x, y, &value);
			if(value < threshold)
				value = 0;
			else
				value = 255;
			FreeImage_SetPixelIndex(chg, x, y, &value);
		}
	}
}

/* 图像标定
 * 让灰度值分布均匀
伽马校正后，一定需要再次标定。
所以，伽马校正时，先存储归一化的float值；
校正结束后，再恢复到 [0, 255]
这样在伽马值大于1时，也很方便。否则很容易溢出
value = 255 * (v - min) / max
*/
static void scale_color(FIBITMAP* dib, void* data, float min, float max)
{
	//val = 255 * (val - min) / max
	int x, y, index;
	BYTE value;
	float v;

	int w, h;
	w = FreeImage_GetWidth(dib);
	h = FreeImage_GetHeight(dib);

	float *ptr = (float*)data;
	for (y = 0; y < h; y++) {
		for (x = 0; x < w; x++) {
			index = get_window_index(y, x, w);
			v = ptr[index];
			value = 255 * (v - min) / max;
			FreeImage_SetPixelIndex(dib, x, y, &value);
		}
	}
}

static void sinit_extreme(float* min, float* max)
{
	*min = 1.0f;
	*max = 0.0f;
}

static void sget_extreme(float value, float* min, float* max)
{
	if (value < *min)
		*min = value;
	if (value > *max)
		*max = value;
}

/*
 * 伽马校正，data用于存放中间结果
 * 归一化后，先写入data, 再写入 chg
 * */
void img_gamma(FIBITMAP* ori, FIBITMAP* chg, void* data, float gamma)
{
	int x, y, index;
	BYTE value;

	int w, h;
	w = FreeImage_GetWidth(ori);
	h = FreeImage_GetHeight(ori);

	float min, max;
	float v;
	float *ptr = (float*)data;
	sinit_extreme(&min, &max);
	for (y = 0; y < h; y++) {
		for (x = 0; x < w; x++) {
			FreeImage_GetPixelIndex(ori, x, y, &value);
			v = (float)value / 255; //归一化
			v = powf(v, gamma); //伽马校正，归一化后的值，幂次值方便计算
			sget_extreme(v, &min, &max); //更新极值
			index = get_window_index(y, x, w);
			ptr[index] = v;
		}
	}
	scale_color(chg, data, min, max);
}

/*
 * 计算所有灰度值之和
 * */
unsigned int histo_get_total(unsigned int* histo)
{
	unsigned int total = 0;
	int i;
	for (i = 0; i < 256; i++)
		total += histo[i];
	return total;
}

/* 计算均值 */
unsigned char histo_get_image_avg(FIBITMAP* dib, unsigned int* histo)
{
	unsigned int total, avg;
	int w, h, x, y;
	BYTE value;
	w = FreeImage_GetWidth(dib);
	h = FreeImage_GetHeight(dib);
	
	total = 0;
	for (y = 0; y < h; y++) {
		for (x = 0; x < w; x++) {
			FreeImage_GetPixelIndex(dib, x, y, &value);
			total += value;
		}
	}
	avg = total / w / h;
	return (unsigned char)avg;
}

//计算整幅图灰度的标准差
unsigned int histo_get_image_sigma(FIBITMAP* dib, unsigned char avg)
{
	int w, h, x, y;
	w = FreeImage_GetWidth(dib);
	h = FreeImage_GetHeight(dib);

	unsigned int sigma;
	unsigned int total = 0;
	BYTE value;
	for (y = 0; y < h; y++) {
		for (x = 0; x < w; x++) {
			FreeImage_GetPixelIndex(dib, x, y, &value);
			total += (value - avg) * (value - avg);
		}
	}
	sigma = total / w / h;
	sigma = (unsigned int)sqrt(sigma);
	return sigma;
}

//unsigned char sk[256]
void get_CDF(unsigned char* sk, float* prk)
{
	int i;
	float temp = 0.0f;
	for (i = 0; i < 256; i++) {
		temp += prk[i];
		sk[i] = 255.0f * temp;
	}
}

//calc sk
void get_sk_value(unsigned int* histo, unsigned char* sk, int wh)
{
	int i;

	//prk
	float prk[256];
	for (i = 0; i < 256; i++) {
		if (histo[i])
			prk[i] = (float)histo[i] / wh;
		else
			prk[i] = 0.0f;
	}

	get_CDF(sk, prk);
}

/* 在gzq[i] 中找到 sk中的值，并记录下来 */
void histo_get_new_mapping(unsigned char* gzq, unsigned char* sk)
{
	int i, j;
	unsigned char value;
	for( i = 0; i < 256; i++) {
		value = sk[i];
		for(j = 0; j < 256; j++) {
			if(gzq[j] >= value) {
				sk[i] = j;
				break;//found it
			}
		}
	}
}

void histo_apply_mapping(FIBITMAP* ori, FIBITMAP* chg, unsigned char* sk_data)
{
	int x, y;
	int w, h;
	BYTE value, new_value;
	w = FreeImage_GetWidth(ori);
	h = FreeImage_GetHeight(ori);

	for (y = 0; y < h; y++) {
		for (x = 0; x < w; x++) {
			FreeImage_GetPixelIndex(ori, x, y, &value);
			new_value = sk_data[value];
			FreeImage_SetPixelIndex(chg, x, y, &new_value);
		}
	}
}

/*直方图均衡 
 * 变换函数存储在 sk_data[256] 中*/
void image_histogram_equa(FIBITMAP* ori, FIBITMAP* chg, unsigned char* sk_data)
{
	int w, h, wh;
	unsigned char sk[256];

	w = FreeImage_GetWidth(ori);
	h = FreeImage_GetHeight(ori);
	wh = w * h;

	unsigned int ori_his[256];
	get_histogram(ori, ori_his);
	get_sk_value(ori_his, sk, wh);
	if(sk_data)
		memcpy(sk_data, sk, 256);
	//use new map
	histo_apply_mapping(ori, chg, sk);
}

/*
* 二次高斯拟合
*/
void get_match_histogram(unsigned int* chg_histo)
{
	//求取指定直方图的CDF
	int i;
    for(i = 0; i < 256; i++){
		chg_histo[i] = twoModeGauss(i);
    }
	//print_histo_data(chg_histo);
}

/*
* 统计 (x, y) 为中心点的局域直方图
* e.g 3 * 3 方框内的直方图
* 
*/
void calc_local_his(FIBITMAP* dib, int x, int y, int type, unsigned int* histo)
{
	int i, j;
	int r = type / 2;
	BYTE val;
	for (i = y - r; i <= y + r; i++) {
		for (j = x - r; j <= x + r; j++) {
			FreeImage_GetPixelIndex(dib, j, i, &val);
			histo[val]++;
		}
	}
}

/*
 * 直方图均衡化将该邻域
 * type = 3, 3 * 3 邻域
 * type = 5, 5 * 5 邻域
 */
void image_histogram_local(FIBITMAP* ori, FIBITMAP* chg, int type)
{
	int w, h;
	w = FreeImage_GetWidth(ori);
	h = FreeImage_GetHeight(ori);
	
	unsigned int histo[256]; //hitogram
	unsigned char sk_data[256]; //sk_data
	unsigned char new_value, sub;
	int r = type / 2;
	int area = type * type;
	
	//第3步：移动邻域，扫描整个图像
	int y, x;
	for(y = r; y < (h - r); y++){
		for(x = r; x < (w - r); x++) {
			memset(histo, 0, sizeof(histo));
			calc_local_his(ori, x, y, type, histo);
			get_sk_value(histo, sk_data, area);
			FreeImage_GetPixelIndex(ori, x, y, &sub);
			new_value = sk_data[sub];
			FreeImage_SetPixelIndex(chg, x, y, &new_value);
		} //end for x
	}//end for y
}

/*统计直方图
 * type: 3 * 3, 5 * 5
 * avg, sigma: 全局均值，标准差
 *
 ** if(k0 * avg >= lavg && sigma_l>=k1*sigma_g && sigma_l <= k2*sigma_g)
 *		新的灰度值 = e * 原图的灰度值
 * avg: 全局图像均值
 * lavg: 局部均值
 * sigma_g: 图像的标准差
 * sigma_l: 局部标准差
 * */
void image_histogram_stats(FIBITMAP* ori, FIBITMAP* chg, int type,\
	unsigned char avg, unsigned int sigma,\
	float e, float k0, float k1, float k2)
{
	float lavg;
	float lsigma;

	int y, x;
	int w, h;
	int p, q;

	int size = type * type;
	unsigned char *histo = malloc(size);

	w = FreeImage_GetWidth(ori);
	h = FreeImage_GetHeight(ori);
	int dis = type / 2;
	int k;
	BYTE val;
	for(y = dis; y < (h - dis); y++){
		for(x = dis; x < (w - dis); x++) {
			memset(histo, 0, size);
			k = 0;
			for(p = x - dis; p <= x + dis; p++) {
				for(q = y - dis; q <= y + dis; q++) {
					FreeImage_GetPixelIndex(ori, p, q, &histo[k]);
					k++;
				}
			}
			lavg = data_get_avg(histo, size);
			lsigma = data_get_sigma2(histo, size, lavg);
			lsigma = sqrtf(lsigma);
			if(lavg <= k0 * avg && lsigma <= k2 * sigma &&\
				lsigma >= k1 * sigma) {
				//update
				FreeImage_GetPixelIndex(ori, x, y, &val);
				val = e * val;
				FreeImage_SetPixelIndex(chg, x, y, &val);
			}//end if
		} //end for x
	}

	free(histo);
}


